# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2020/12/9 14:16
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : norm.py
@Version     : Version 1.0.0
@Description : TODO 自己实现计算向量范数，矩阵范数，矩阵条件数
@Created By  : PyCharm
"""
import numpy as np
from scipy.linalg import eigvals, det, inv


def vector_norm(vector: np.ndarray, p=None):
    """
    计算向量的p-范数
    :param vector: 实向量或者复向量
    :param p: 指定类型的范数，默认是oo范数
    :return: 指定向量范数和p值
    """
    if p is None:
        return abs(vector).max(), p
    elif p >= 1:
        return np.power(np.sum(np.power(abs(vector), p)), 1 / p), p
    else:
        raise Exception("error,p must be an integer , greater than  or equal to 1")


def matrix_norm(square_matrix: np.ndarray, v=None):
    """
    计算矩阵范数，其中v=1,2,oo
    :param square_matrix: 实数方阵或者复数方阵
    :param v: 范数类型，默认是弗罗贝尼乌斯范数
    :return: 指定类型范数
    """
    if square_matrix.shape[0] == square_matrix.shape[1]:
        if v is None:
            # 弗罗贝尼乌斯范数
            return np.power(np.sum(np.power(square_matrix, 2)), 1 / 2)
        elif v == 0:
            # 矩阵的列范数
            return np.max(np.sum(abs(square_matrix), v))
        elif v == 1:
            # 矩阵的行范数
            return np.max(np.sum(abs(square_matrix), v))
        elif v == 2:
            # 矩阵的2-范数
            return np.sqrt(eigvals(square_matrix.T @ square_matrix).max()).real
        else:
            raise Exception("sorry,currently,v must be None,0,1 or2")
    else:
        raise Exception("error,input matrix must be a square matrix.")


def condition_number(square_matrix: np.ndarray, v=None):
    """
    计算矩阵条件数，其中v=oo,1,2
    :param square_matrix:实数方阵或者复数方阵
    :param v:范数类型，默认是oo范数
    :return:指定范数下的矩阵条件数
    """
    if square_matrix.shape[0] == square_matrix.shape[1]:
        if det(square_matrix) != 0:
            if v is None:
                return matrix_norm(square_matrix, v=1)*matrix_norm(inv(square_matrix), v=1)
            elif v == 1:
                return matrix_norm(square_matrix, v=0)*matrix_norm(inv(square_matrix), v=0)
            elif v == 2:
                eigenvalue_fraction = matrix_norm(square_matrix, v=2) / np.sqrt(
                    eigvals(square_matrix.T @ square_matrix).min()).real
                return np.sqrt(eigenvalue_fraction)
            else:
                raise Exception("sorry,currently,v must be None,1 or2")
        else:
            raise Exception("error,input matrix must be a non-singular matrix")
    else:
        raise Exception("error,input matrix must be a square matrix.")


if __name__ == '__main__':
    real_array = np.array([1, -2, 3]).reshape(3, 1)
    complex_array = np.array([[1 + 2j, -3 - 5j, 4 + 1j],
                              [7 + 8j, 5.2, 6.7],
                              [3 - 2j, 1 - 1j, 6.3 + 4j]],
                             dtype=np.complex128)
    # 计算矩阵范数，来源详见李庆扬数值分析第5版P166
    """
    a = np.array([[1, -2], [-3, 4]])
    row_norm=matrix_norm(a,v=1)
    column_norm=matrix_norm(a,v=0)
    matrix_second_norm=matrix_norm(a,v=2)
    frobenius_norm=matrix_norm(a)
    print(column_norm,row_norm,frobenius_norm,matrix_second_norm,sep='\n')
    """
